Nonlinear screening and percolative transition in a two-dimensional electron liquid 
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A variational method is proposed for calculating the percolation threshold, the real-space struc- 
ture, and the ground-state energy of a disordered two-dimensional electron liquid. Its high accuracy 
is verified against exact asymptotics and prior numerical results. The inverse thermodynamical den- 
sity of states is shown to have a strongly asymmetric minimum at a density that is approximately the 
triple of the percolation threshold. This implies that the experimentally observed metal-insulator 
transition takes place well before the percolation point is reached. 



The discovery that a two-dimensional (2D) electron 
Hquid can be a metal at moderate electron density rie 
and an insulator at small rig was a major surprise that 
questioned our fundamental understanding of the role of 
disorder in such systems. Today, almost a decade later, 
it remains a subject of an intense debated One impor- 
tant reason why the conventional theory fails could be its 
flawed basic premise of the "good" metal, i.e., a uniform 
electron liquid slightly perturbed by impurities and de- 
fects. Indeed, modern nanoscale imaging techniques2i2i^ 
unambiguously showed that low-density 2D electron sys- 
tems are strongly inhomogeneous, "bad" metals, where 
effects of disorder are nonperturbatively strong. In par- 
ticular, depletion regions (DR), i.e., regions where n(r) is 
effectively zero, exist. They appear when rig is too small 
to adequately compensate fluctuating charge density of 
randomly positioned impurities. As Tig is reduced in the 
experiment, e.g., in order to approach the vicinity of the 
metal-insulator transition (MIT), the DRs are expected 
to grow in size and concentration and eventually merge 
below some percolation threshold Ue = rip. An important 
and controversial issue is whether or not this percolation 
transition plays any role in the observed MITi To resolve 
it one needs to have a theory that is able to calculate 
Up and that can describe the inhomogeneous structure 
of the 2D metal at Tie ~ rip. A great progress in this 
direction has been achieved by Efros, Pikus, and Bur- 
nett whose paper is intellectually tied to earlier work 
on nonlinear screening by Efros and Shklovskii.^ Still, 
analytical results remained scanty and numerical simula- 
tionsSi^ were the only known way to quantitatively study 
the rif, ~ Up regime. These simulations are very time 
consuming and redoing them in order to get any infor- 
mation beyond what is published^'* or to study novel 
experimental setups seems impracticable. Below I will 
show that a variational approach to the problem can be 
a viable alternative. Comparing it with the available nu- 
merical results for a typical model of the experimental 
geometry (Fig. , I establish that it correctly predicts 
the value of Up and accurately reproduces the energetics 
of the ground state, in particular, the density dependence 
of the electrochemical potential and of the inverse ther- 
modynamical density of states (ITDOS), = d^/due, 
over a broad range of rig. The most striking feature of 
the resultant function x~^("-e) is a strongly asymmetric 
minimum, which is observed in real experimentsAiS*iiii^ 



FIG. 1; The geometry of the theoretical model. The 2D 
layer of interest is sandwiched between the top and the bot- 
tom metallic gates. The dopants reside in the plane with the 
dashed border. The depletion regions (shown as holes in the 
probed layer) enhance the penetrating electric field Ep. 



I will elaborate on the origin of this feature and show 
that it occurs at rig « Sup largely independently of the 
parameters of the system. Recently, this minimum at- 
tracted much interest when Dultz et a/.^° reported that 
in some samples it virtually coincides with the appar- 
ent MIT. The proposed theory indicates that, at least 
for these samples, any connection between the MIT and 
the percolation of the DRs can be ruled out. Thus, the 
explanation of the MIT lies elsewhere. 

The Hamiltonian of the model is adopted from Efros, 
Pikus, and Burnett (EPB)^ (see also Fig.Ul), 



H = I d^r 



1, 



[n(r) -ne]$(r) + i/o(n) 



where <I>(r) is the electrostatic potential. 
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K is the dielectric constant, and Ho(n) is the energy den- 
sity of the uniform liquid of density n, 



(3) 



At low densities, ho{n) ^ 1 is a slow functioni^ of 
n. The negative sign in Eq. Q reflects the preva- 
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lence of the exchange-correlation energy over the ki- 
netic one in this regime. The potential created by 
the random concentration of dopants nrf(r) can be ex- 
pressed in terms of the effective in-plane background 
charge o'(q) = '^^(q) exp(— gs) (the tildes denote the 
Fourier transforms). With these definition, $(r) coin- 
cides with the potential created by the charge density 
cr(r) = n(r) — Ue — o'(r). I will assume that the dopants 
have the average concentration (n^) = 3> and are 
uncorrelated in space. In this case the one- and two-point 
distribution functions of <j{v) have the Gaussian form. 



A (a) 
1 



(2^i^o)-i/2exp(-aV2i^o), 



exp 



2Kr'(7cr' — 



l2\-\ 



(4) 
(5) 



where cr' = (T(r r') and Kr is given by 



Kr 



(CT(r)CT(O)) = n,s/Ti{r^ + 45^)3/2. 



(6) 



The characteristic scales in the problem are as follows.^ 
The typical amplitude of fluctuations in a is yjulj see 
Eqs. (@J and ((HJ. Their characteristic spatial scale is the 
spacer width s [cf. Eqs. £^nd Fig. Q]. In the cases 

studied below, rig ^ Up and rig 3> Up, s exceeds the 

— 1/2 

average interelectron separation ao = Ue ■ As usual 
in Coulomb problems, the energy is dominated by the 
longest scales, in this case s. The fluctuations Ho{n) — 
Ho{ne) of the the local energy density come from the 
interactions on the much shorter scale of oq <SC s and 
can be treated as a perturbationi& For the purpose of 
calculating the ground-state density profile n(r) I neglect 
Hq. Once such a ground state is known, I correct the 
total energy H by adding to it Hq averaged over n(r). 

To find n(r) one needs to solve the electrostatic prob- 
lem with the following dual boundary conditions: if 
7T.(r) > 0, then $(r) = const; otherwise, if rt = 
(DR), thenSiLii ^ > fi. I start with the analysis of the 
large-rie case, which clarifies why x~^{''^e) dependence is 
nonmonotonic and which provides a formula for the the 
density rim where has the minimum. 

In the limit Ue 3> yjulj the asymptotically exact 
treatment is possible because the DRs appear only in 
rare places where cr(r) dips below — rig. (In practice, this 

limit is realized when Uf. > K^^ « ^.I^^Jnlj s). The cor- 
responding electrostatic problem is analogous to that of 
the metallic sheet perforated by small holes, see Fig. ^ 
The most elegant way to derive x~^ iii this regime is 
to calculate the fraction EpjE^ of the electric field that 
reaches the bottom layer in the geometry of Fig. ^ In- 
deed, Ep/Eq is nonzero only if the probe layer is not a 
perfect metal, x~^ 7^ 0- I^i the simplest case, where dis- 
tances si and S2 are large, the following formula holds: 
X'^ = ATT{e^S2/K){dEp/dEo) (cf. Ref.l). This is es- 
sentially the formula used to deduce x~^ in the exper- 
imen^MiiSiiiii^) . It is immediately obvious that holes in 
the metallic sheet enhance the penetrated field Ep. For 
example, the field leaking through a round hole of radius 



a is the field of a dipolei^ p = (a3/37r)Eo. If there is a fi- 
nite but small concentration Nf^ of such holes, their fields 
are additive, leading to x^^ = {Sne'^ /3K,)Nha^ ■ From 
here the exact large-ng asymptotics of x~^ and /i can be 
obtained by substituting the proper and averaging 
over the distribution of a. This can be done by noting 
that these holes appear around the minima of (T(r) whose 
statistics is fixed by Eqs. It is easy to see that the 

most probable holes are nearly perfect circles with radii^ 
a ~ s^/Kq/ui,. The charge distribution around a single 
hole at distances r > a is given by the formula 



n{r) 



/ r*^ 2 \ a 
— — 1 — — + - I arccos - 



, (7) 



where the hole is assumed to be centered at r = and 
= — 3[ne + (T{0)]/axx- Equation {Tj) can be obtained, 
e.g., by generalizing the textbook solution^^ for the hole 
in the metallic sheeti^ and is also the limiting form of 
Eq. (11) in Ref. [H Previouslj^ Eq. was used for 
study of quantum dots in Ref. Il8l 

In the current problem the main factor that determines 
the net contribution Xdr of the depletion holes to x^^ is 
their exponentially small concentration, proportional to 
Pi(— Tig) cx exp{—nl/2Ko). The final result. 



{i^/2/STT){e^n,/Ksnl) exp(-47rs2n^/n,), (8) 



agrees with that of EPB but has no numerical coefficients 
left undetermined. To finalize the calculation, one needs 
to augment Xdr by the local term Xo^ = cP{Ho{n)) /dnl. 
In the present case, fluctuations around the average den- 
sity are small. Hence, {Ho{n)) ~ iJo("-e) and 

Xo^iue) ~ -{e^/K)hi{ne)/y/n^, ne:>y/ru/s, (9) 



where /ii(ne) = (3/4)/io(rie) + S/iprie + h'^nl ~ 1. 

Formula © implies a sharp upturn of the ITDOS as 
Tie ^ caused by the exponential growth of the DRs. 
At the boundary of its validity, rie ~ ^/tth/s, Eq. ^ 
gives Xdr ~ b'^s/k, which is large and positive. On the 
other hand, the local term Xo^ i^'i- ®] shows a weak 
dependence on n^, remaining small and negative. Com- 
bined, they produce a strongly asymmetric minimum in 
X"^("-e) = Xdr + Xo^ at the density 



4^ S [TThj ' 



(10) 



The log- factor in Eq. (|10l) is rather insensitive to n^s^ 
and hi. For n^s^ — hi = 1, one gets Um ~ 0.38y/ni/ s. 
In comparison, (see Ref. 6 and below) the percolation 
threshold is 



Up « 0.12^/n7/ s. 



(11) 



so that 



3n„. In accord with EPB's heuristic ar- 



gument, at such density a very small area fraction is 
depleted (see inset in Fig. |2t), and so the use of the 
asymptotic formula (|SJ| is justified. 
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FIG. 2: (a) Electrochemical potential vs. electron density for 
the electrostatic problem, Hq = 0, according to the present 
theory (solid line) and EPB's numerical simulations (dots). 
Inset: Fraction of the depleted area vs. density, (b) ""S. 
density according to the numerical simulations of Shi and Xie^ 
(squares), the present theory (solid line), and for the uniform 
electron liquid (dashed line). 



Let us now proceed to the case ~ rip. Again, I start 
with the electrostatic problem (iJo = 0). One expects 
DRs to be abundant and irregularly shaped. For a given 
rie, the ground-state n(r) is some nonlocal functional of 
a and there is no hope to find it exactly. What I wish to 
report here is that a variational solution sought within 
the class of purely local functionals, n(r) = n[f7(r)] re- 
markably accurately reproduces the ^{rie) and x~^("-e) 
dependencies found in numerical work.^'- More interest- 
ingly, it predicts the correct value of Up [Eq. ((TT|l ] . 

The system of equation that defines such a variational 
state n{a) follows from the fact that for a Gaussian ran- 
dom function the averages over the total area of the 
system and over the distribution function are the same. 
This yields [cf. Eqs. 



G(ct, a') - 



dcrda'p{cr)G{a,a')p{a'), (12) 
d\'V{r')[P2{a,a') - Pi(a)Pi(a')], (13) 



energy Hy needs to be minimized with respect to all func- 
tions n(a) that obey the constraints n{a) > and 



d<j Pi{a)n{ij) 



(14) 



The latter ensures that the average density is equal to 
rie. Introducing the Lagrange multiplier /i^ (variational 
estimate of the electrochemical potential), one obtains 
that Hy is minimized if, for all a > f, p{a') satisfies 



d(j'G{a,a')p{a') = py{ne)Pi{(y) 



(15) 



Here / is such that n(/) = 0. [I found that n{a) is al- 
ways a monotonically increasing function, so that n > Q 
corresponds to cr > /]. The kernel G{a,a') [Eq. H13|l ] is 
log-divergent, G cx — Injcr — a'| at a a' and decays 
exponentially at large a and cr'. There is a certain anal- 
ogy between Eq. H15|) and the integral equations of ID 
electrostatics, which also have log-divergent kcrnelsii^ii^ 
This analogy entails that n{(j) ^ ^/c — f at cr close to /. 
Note that cr(r) — / is proportional to the distance in the 
real space between the given point r and the boundary 
of the nearby DR. Thus, the variational principle renders 
correctly the square-root singularity in n(r) at the edge 
of the DRs [cf., e.g., Eq {Tj)]. I was not able to establish 
the analytical form of the solution beyond this property 
and resorted to finding n{a) numerically. To do so the 
integral in Eq. (|15|l was converted into a discrete sum 
over 101 points on the interval |cr| < LSy^/s the re- 
sultant system of 101 linear equations was solved on the 
computer. The solution can be approximated by a simple 
analytical ansatz 

naicT) - [{n, + af - [n, + !f\^lH{cy - /), (16) 

where d{z) is the step-function. For example, the energy 
Hy is nearly the same whether it is calculated using Ua 
or using the actual solution of Eq. (|15|) . So, in principle, 
Eq. (|16l) obviates the need to solve Eq. H15|) . The only 
equation that needs to be solved is Eq. p4|l for /. 

At this point one can compare the predictions of the 
variational method for /z„(ne) with EPB's numerical re- 
sults that were also obtained for the Hq = problem. As 
Fig. EK illustrates, they are in a good agreementpi^ 

To test the theory further I compare it next with the 
numerical data of Shi and Xiei^ To this end, one needs 
to calculate including the effect of a finite Hq. The 
first step is to take the derivative Xdr ~ dpy/due, which 
is easily done numerically. An accurate fit to the result 
is provided by the interpolation formula 



9 



?,^/2 0.30 + 7/ 



87r?7 0.036 + 0.12?7 -f?7' 



■exp(-47r77^), (17) 



where /o(cr) = n{a) — — a and V{r') 



'-/kt'. The 



where 77 = n^s/ ^Jnl. This particular form is devised to 
match Eq. ||SJ) at large Ue and is consistent up to log- 
corrections with the behavior expectedi^ at very small 
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rie. The next step is to evaluate the quantity 

^^'^d^^^'^^'di^ daPi(a)i7oKa)], (18) 

which is also easily done on the computer. The total 
ITDOS, x~^{ne) = XDR + Xo^ calculated from Eqs. ifTTl - 
I18|l and the theoretical value of n,„ ~ 0.7 x 10~^a~^^ [per 
Eq. H10() ] compare very well with the simulations, see 
Fig-Eb- The parameters I used are — 6.25 x 10~'*a^^ 
and s — lOag, where is the effective Bohr radiusi^ In 
agreement with the exact results presented above, Xo^ 
X-DR ^ ^^'^ so ^^^^ upturn of x~^ at low Ue is 

driven the growth of DRs. 

Within the variational method the boundaries of the 
DRs coincide with the level lines a{r) = / of the zero- 
mean random function a. Consequently, the DRs perco- 
late at / > 0. Solving Eq. UHJl with / = 0, one arrives 
at Eq. which is in excellent agreement with EPB's 
result Up « 0.11y^/s. One important quantity not re- 
ported in the published numerical worka^ is the area 
fraction of the DR. Within the variational method, it is 
equal to erfc(— //V2Xo)/2, where erfc(z) is the comple- 
mentary error function. It is exactly 1/2 at rig = rip and 
increases as — > as shown in Fig.[2ti (inset). 

Let us now compare our results with the experimen- 
tal data. Taking a rough number = 3 x lO^^cm"^ 
and a typical spacer width s = 40 nm in Eq. (|10|l . 
one gets Um ~ 5.2 x 10^*^ cm^^, in agreement with ob- 
served valuesiiii The estimate for the percolation point is 



Hp « nm/3 « 1.7 X 10^°cm~^. Despite some small uncer- 
tainty in the last number (due to the uncertainties in rii 
and hi), rip and Um differ substantially, and so they can 
be easily distinguished experimentaly. Therefore, the ob- 
served apparent MIT^" at n — n,n has nothing to do with 
the percolation of the DRs and moreover with breaking 
of the electron liquid into droplets (the latter occurs at 
rie <C Up). From Fig. O^, one can estimate that the DR 
occupy mere 6% of the total area at n = rim- 

One qualitative prediction that follows from Eqs. © 
and H17|l is that the upturn of x~^{^e) should be sharper 
in samples with larger s, which seems to be the case if 
the dat a of Ref. M (s — 14 nm) are compared with those 
of Refs. Iinllllia f.s < 4nm). Detailed fits are left for 
future. 

I conclude with mentioning some other theoretical 
work on the subject. It was suggested^^ that near the 
MIT, function x~^('^e) may contain both regular and sin- 
gular parts. My theory can be considered the calculation 
of the former. Its detailed comparison with experiment 
may furnish an estimate of the putative singular term. 
The effect of disorder on x~^ was also studied in Refs.l2^ 
and 13 but the electron density inhomogeneity was not 
accounted for. Finally, a nonmonotonic x~^('^e) depen- 
dence was found in a model^ without disorder. 
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